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Abstract 

This  paper  presents  analytical  correlations,  which  predict  the  polarization  performance  and  thermal  effects  in  an  intermediate  temperature  proton 
exchange  membrane  fuel  cell  (PEMFC).  Such  correlations  are  useful  for  engineers  and  designers  of  fuel  cells  to  expedite  calculations  without 
depending  on  complex  computational  models.  Analytical  results  compare  well  with  published  experimental  polarization  data.  They  also  indicate 
the  temperature  variations  and  dominant  modes  of  heat  transfer  in  a  unit  cell. 
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1.  Introduction 

Low  temperature  (<90  °C)  proton  exchange  membrane  fuel 
cells  (PEMFCs)  have  long  been  considered  as  candidates  for 
future  energy  systems.  However,  numerous  debilitating  factors 
-  such  as  sluggish  electrode  kinetics,  susceptibility  of  platinum 
(Pt)  catalysts  to  CO  poisoning,  and  water  management  prob¬ 
lems  associated  with  Nafion®  membranes  -  have  hindered  their 
progress. 

Much  research  has  been  conducted  on  new  proton  conduct¬ 
ing  membranes  that  are  capable  of  operating  at  intermediate 
temperatures  (150-200  °C)  [1].  At  these  temperatures  CO  poi¬ 
soning  becomes  less  prominent,  electrode  kinetics  are  faster, 
and  water  would  exist  primarily  in  the  vapor  phase  precluding 
problems  associated  with  water  management  and  mass  transport 
limitations. 

Although  modeling  of  PEMFCs  has  received  much  attention 
over  the  past  15  years  [2],  most  of  the  work  is  based  on  PEM¬ 
FCs  using  a  Nafion®  membrane.  In  a  previous  paper  [3],  we 
developed  a  mathematical  model  for  a  PEMFC  using  a  polyben¬ 
zimidazole  (PBI)  membrane.  Mathematical  models,  however, 
tend  to  be  very  complex  and  esoteric,  and  do  not  parametrically 
illustrate  the  dependence  of  the  final  solution  on  factors  affecting 
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it.  They  provide  detailed  and  accurate  calculations  for  complex 
designs,  however,  analytical  correlations,  even  if  approximate, 
enable  quick  calculations  for  most  engineering  purposes. 

This  paper  presents  correlations  based  on  analytical  solutions 
for  the  model  of  the  PBI  PEMFC  developed  in  ref.  [3].  The  ana¬ 
lytical  results  are  compared  to  the  experimental  results  published 
by  Wang  et  al.  [4]. 

2.  Mathematical  formulation 

The  domain  of  interest  is  the  membrane  electrode  assembly 
(ME A).  For  the  purpose  of  a  ID  study,  the  gas  flow  channels  are 
not  considered;  instead  boundary  conditions  are  derived  at  the 
gas  channel/gas  diffuser  interfaces.  Fig.  1  shows  the  schematic  of 
the  membrane  electrode  assembly  (MEA)  and  the  dimensionless 
coordinates  used. 

Only  steady  state  operation  is  considered,  and  at  150  °C, 
water  is  expected  to  exist  only  in  the  vapor  phase,  and  all  electro¬ 
chemical  reactions  occur  in  the  gaseous  phase.  The  gas  diffuser 
regions  are  modeled  as  macro-homogenous  porous  media.  The 
catalyst  is  treated  as  an  interface,  where  species  are  consumed 
and  produced. 

2.7.  Equations 

The  phenomenological  equations,  based  on  the 
Stefan-Maxwell  equations,  Darcy’s  law,  conservation  laws, 
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Nomenclature 

A  ME  A  cross  sectional  area 

c  specific  heat  capacity 

D  gas  pair  diffusivity 

E  potential 

F  Faraday  constant,  96,487  C  mol-1 

I  current  density 

k  thermal  conductivity 

/  thickness 

m  mass  fraction 

M  molecular  mass 

N  mass  flux 

P  pressure 

R  universal  gas  constant,  8.3143  J  mol-1  K-1 
S  source,  entropy 

T  temperature 

u  mass  averaged  velocity 

V  volumetric  flow  rate,  potential 

Greek  letters 

a  charge  transfer  co-efficient 

/3  porous  media  gas  permeability 

s  gas  porosity 

fi  dynamic  viscosity 

p  density 

a  conductivity 
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Fig.  1.  Schematic  of  the  MEA  with  dimensionless  coordinates. 


The  entropy  change  term  for  the  overall  fuel  cell  reaction  with 
water  in  the  vapor  form  is  correlated  using  thermodynamics 
tables.  The  following  relation  is  valid  for  400  K  <  T<  500  K 

—  =  -18.449 -0.01283TJ  mol-1  K_1  (3) 

n 

The  activation  overpotential  at  the  cathode  catalyst  layer  is 
calculated  from  the  Tafel  relation 


RT 

2aF 


In 


and  ideal  gas  laws,  are  shown  in  Table  1,  as  derived  in 
ref.  [3]. 

Once  the  partial  pressures  of  hydrogen,  oxygen  and  vapor 
at  their  respective  catalyst  layers  are  determined,  the  reversible 
cell  potential  at  a  given  temperature  can  be  calculated  from  Eq. 
(1).  The  reversible  cell  potential  at  Tref  (=298  K)  is  1.185  V  for 
water  existing  in  the  vapor  phase 


Tref)  + 


2  RT 
nF 


In 


m2  x  p0o52 


P  w 


The  anode  overpotential  is  neglected  since  at  elevated  tem¬ 
peratures,  the  electrode  kinetics  associated  with  hydrogen  disso¬ 
ciation  is  expected  to  be  much  faster  than  for  oxygen  reduction. 
The  total  ohmic  overpotential  is  given  by 


T^ohm  —  7 


The  operating  cell  potential  is  given  by  the  reversible  poten¬ 
tial  (concentration  effects  included)  less  the  activation  and  ohmic 
overpotentials 


where 


rrii/Mi 

Pi  ~  PJ2mi/Mi 


T^cell  —  Fact  7s0h  m  (6) 

The  energy  equation  over  the  entire  unified  domain  is  given 
by  Eq.  (7),  which  accounts  for  convection,  conduction  and  heat 


Table  1 

Summary  of  governing  equations 


Anode 


Fluxes 

Species  conservation 

Ideal  gas  relation 

Mass  conservation 
Darcy’s  law 
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generation.  In  the  membrane,  since  there  is  no  fluid  flow,  the 
convection  term  vanishes 


d  (cT)  ,  d  2T 

PU — ; =  &eff“r^r  + 


d  z 


d  z2 


(7) 


The  heat  generation  is  due  to  ohmic  heating  and  heat  of  reaction 


—  ^ohmic  T  ^reaction 


5ohmic  — 


a 


(8) 

(9) 


The  heat  of  reaction  is  considered  to  be  concentrated  in  the 
cathode  catalyst  layer.  It  is  due  to  the  entropy  of  reaction  and 
the  irreversibility  associated  with  activation  overpotential 


’reaction 


=  / 


(10) 


2.2.  Boundary  conditions 


2.3.  Analytical  solutions 


It  can  be  shown  that  the  solution  of  the  general  ordinary  dif¬ 
ferential  equation 


dy(v) 

=  A(x)(B  +  Cy)(D  +  Ey) 
ax 

with  boundary  condition  y(x\)  =  y\,  and  DC ^ BE ,  is 

DH-B 

y  = - 

7  C-  EH 


(15) 


(16) 


where 


H  = 


B  +  Cy  i 
D  +  Eyi 


exp  [(DC  —  BE)A(x  —  jq)] 


A(x-x i)=  f  A(x) dx 


(17) 

(18) 


The  Stefan-Maxwell  equation  can  be  combined  with  the  ideal 
gas  relationship  to  obtain  an  ODE  similar  to  Eq.  (15) 


At  the  interface  of  the  gas  flow  channel  and  gas  diffusers,  the 
mass  fraction  boundary  conditions  can  be  determined  from  a 
material  balance  between  gas  supply  and  reaction  requirements 


at  f  =  0, 
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at  f  =  3, 
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(12) 


at  f  =  3, 

mN2  =  — o 

(mN2,IN 


_ mN2,IN _ 

Axs|A^o2  I)  T  (^n2,in)  (^w,IN  A-xs l^w I) 


(13) 


Since  each  gas  occupies  the  same  volume,  the  inlet  mass  flow 
rate  for  each  species  can  be  calculated  from  Eq.  (14) 


m 


o  _ 

z,IN  — 


RT 


(14) 


The  temperature  at  which  supply  gases  are  humidified  deter¬ 
mines  the  partial  pressure  of  the  water  vapor.  The  partial  pres¬ 
sures  of  oxygen  and  nitrogen  are  determined  from  the  stoichio¬ 
metric  make  up  of  air  (21%  oxygen,  79%  atmospheric  nitrogen 
by  molar  composition). 

The  boundary  conditions  for  the  temperature  and  pressure 
are  set  to  the  supply  temperatures  and  pressures  at  the  inter¬ 
faces  between  the  gas  flow  channels  and  the  gas  diffusers.  At 
the  interfaces,  continuity  of  all  variables  is  specified. 
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It  is  assumed  that  due  to  the  low  viscosity  of  the  gas  phases,  the 
pressure  drops  will  be  negligible.  The  mass  fraction  of  hydrogen 
at  the  anode  catalyst  layer  becomes 


ac  _  H  (Mh2/Mw) 

Hz  “  1  -  (Mh2/Mw)  +  H 
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(21) 


In  the  cathode  gas  diffuser,  there  are  two  dependent  species 
variables.  For  simplification,  it  will  be  assumed  that  the  mass 
fraction  of  nitrogen  remains  approximately  constant.  Species 
flow  under  the  influence  of  a  total  pressure  gradient  (convection) 
and  a  concentration  gradient  (diffusion).  Since  the  total  pressure 
gradient  is  expected  to  be  negligible,  the  concentration  gradient 
becomes  the  dominant  factor  affecting  species  flux.  Also,  since 
there  is  no  net  flux  of  the  inert  nitrogen  species  in  the  MEA,  its 
concentration  or  partial  pressure  gradient  is  negligible.  Hence, 
the  mass  fraction  is  expected  to  be  approximately  constant. 

The  corresponding  equation  in  the  cathode  gas  diffuser  is 
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the  solution  of  which,  is 

cc  DH  -  B 
m02  = 


C  -  EH 


where 
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(30) 


The  partial  pressures  necessary  to  determine  the  cell  potential 
can  now  be  found 


P  h2  = 
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All  the  required  parameters  are  available  at  this  point,  except 
the  average  temperature  in  each  diffuser  region.  For  a  simplified 
analysis,  the  temperature  variation  (in  K)  could  be  considered 
negligible,  in  which  case,  T  could  be  replaced  by  Tfi  in  Eqs.  (21) 
and  (24). 

It  is  still  desired  though  to  perform  a  detailed  thermal  analysis. 
The  energy  equation  consists  of  conduction,  convection,  and 
generation  terms.  The  generation  term  can  be  broken  up  into 
ohmic  heating  (volumetric)  and  heat  of  reaction  (interfacial). 
Also  of  note  is  that  in  the  membrane,  there  is  no  fluid  flow, 
hence  the  convection  term  vanishes.  In  a  given  domain  with 
BCs,  T(0)  =  Ti  and  T(L)  =  Tfl,  the  analytical  solution  of  the  energy 
equation  is 
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Table  2 

Parameters  in  Eqs.  (41)  and  (42) 


Full  energy  analysis 
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For  conduction  without  convection,  the  corresponding  equations 
become 


(To  -  Ti)  | 


(38) 


dr 
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Note  that  B  has  the  same  numerical  value  in  both  the  anode  and 
cathode  diffusers,  hence  the  subscript  (d)  is  used.  However,  A 
has  different  values  in  the  anode  and  cathode  diffusers  because 
different  gas  mixtures  will  have  different  specific  heat  capacities. 
Hence,  the  subscripts  (a)  and  (c)  are  used. 

These  results  are  applied  to  the  three  domains  of  the  MEA 
with  continuity  of  temperature  at  each  interface,  and  continuous 
heat  flux  at  the  anode  catalyst  interface.  At  the  cathode  catalyst 
interface  where  the  heat  of  reaction  is  generated,  the  following 
energy  balance  applies 


Reaction  =  (40) 

V  d z/m,2  V  dzyCi2 

The  subscripts  (m)  and  (c)  refer  to  the  heat  transfers  through  the 
membrane  and  cathode  regions,  respectively,  while  (2)  refers  to 
the  position,  f  =  2.  Solving,  we  get 


FG  -  DH 
CF  -  DE 

CH  -  EG 
CF  -  DE 


(41) 

(42) 


Table  2  shows  the  expressions  for  the  above  constants  for 
two  cases:  (1)  the  detailed  energy  analysis  including  conduc¬ 
tion,  convection  and  ohmic  heating,  and  (2)  a  simplified  energy 
analysis  including  conduction  and  ohmic  heating,  but  neglecting 
convection. 
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2.4.  Numerical  values 

The  numerical  values  used  in  the  computation  are  the  same 
as  ref.  [3] .  The  results  presented  in  this  paper  attempt  to  simulate 
the  experimental  results  presented  by  Wang  et  al.  [4].  Inlet  gases 
were  supplied  at  150  °C  (423  K)  and  1  atm.  Both  fuel  (hydro¬ 
gen)  and  oxidant  (oxygen  or  air)  were  humidified  at  28  °C.  So 
we  assume  that  the  partial  pressure  of  water  vapor  in  the  inlet 
streams  is  equal  to  the  saturation  vapor  pressure  at  28  °C. 

3.  Results  and  discussion 

Fig.  2  compares  the  theoretical  IV  curve  with  experimental 
data  for  both  oxidants,  air  and  oxygen.  The  analytical  results 
compare  well  with  the  experimental  data.  For  the  given  oper¬ 
ating  conditions,  respective  peak  power  densities  of  1.58  and 
1.14kW  m-2  for  oxygen  and  air  are  predicted. 

Fig.  3  shows  the  maximum  temperature  increase  at  different 
current  densities  using  air  as  the  oxidant.  The  maximum  tem- 


Fig.  2.  IV  curves  for  analytical  and  experimental  results,  (a)  Analytical  results 
for  oxidant  =  oxygen;  (b)  analytical  results  for  oxidant  =  air. 


Fig.  3.  Maximum  temperature  rise  vs.  current  density,  (a)  Conduction,  convec¬ 
tion  and  ohmic  heating;  (b)  conduction  and  ohmic  heating  with  no  convection; 
(c)  conduction  only  without  convection  and  ohmic  heating. 


2657.1 


2534.5 

◄- 


122.6 


1465.0 

◄ - 


1069.5 


2546.8 

- ► 


2669.4 


4011.8  122.6 


Anode  Gas 

Membrane;  Cathode  Gas 

Diffuser 

Diffuser 

0 

1  \  2 

3  G 

Fig.  4.  Thermal  balance  (W  m  2)  at  I  =  5000  A  m 


.-2 


perature  is  expected  to  occur  in  the  cathode  catalyst  layer  where 
the  heat  of  reaction  is  generated.  The  temperature  increase  is 
then  defined  as  T2  —  Tq.  This  maximum  temperature  variation 
increases  exponentially  with  current  density,  and  for  the  range  of 
current  densities  shown,  is  less  than  2  K.  This  value  represents 
the  ID  temperature  variation  in  the  direction  perpendicular  to 
the  ME  A  at  any  given  section.  This  model  does  not  account 
for  temperature  variations  in  the  direction  of  the  gas  channels. 
Depending  on  the  flow  rate,  the  latter  variation  may  be  signifi¬ 
cant. 

It  can  be  seen  from  Fig.  3  that  neglecting  convection  in  the 
energy  analysis  makes  little  difference  in  the  temperature  varia¬ 
tion.  Thus,  convection  is  a  negligible  means  of  heat  transfer  in  the 
MEA  of  PEM  fuel  cells.  Because  the  solid  material  is  thermally 
conductive  and  the  fluid  velocities  are  very  small,  conduction 
dominates  the  heat  transfer  process. 

Fig.  4  shows  the  magnitudes  of  the  heat  generations  and 
conduction  heat  transfers  at  5000  Am-2.  Ohmic  heating  in  the 
membrane  is  in  the  same  order  of  magnitude  as  the  heat  of  reac¬ 
tion,  but  the  ohmic  heating  in  the  diffusers  are  of  a  smaller  order 
of  magnitude.  This  is  due  to  the  high  electrical  conductivity  of 
the  diffuser  material,  compared  to  the  low  ionic  conductivity  of 
the  PBI  membrane. 

4.  Conclusion 

Analytical  correlations  for  intermediate  temperature  PEM- 
FCs  were  presented.  The  correlations  are  able  to  accurately 
predict  the  polarization  effects,  since  results  compare  favorably 
with  published  experimental  data.  In  addition,  a  ID  temperature 
variation  of  less  than  2  K  is  predicted  within  the  cell  for  typical 
operating  current  densities.  The  results  also  suggest  that  con¬ 
vection  is  not  a  significant  mode  of  heat  transfer  in  intermediate 
temperature  PEMFCs. 
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